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Abstract 

Accurate capturing of discontinuities within compressible flow computations is achieved by cou- 
pling a suitable solver with an automatic adaptive mesh algorithm for unstructured triangular 
meshes. The mesh adaptation procedures developed rely on non-hierarchical dynamical local re- 
finement/derefinement techniques, which hence enable structural optimisation as well as geometrical 
optimisation. The methods described are applied for a number of the ICASE test cases are particu- 
larly interesting for unsteady flow simulations. 

INTRODUCTION 

One of the most important advantages of finite element type unstructured grids is the possibility to 
refine/ derefine locally the mesh during the computation. Successive mesh concentration in “critical” 
zones may be performed, without knowing them a priori during the initial mesh creation, as well as 
coarsening in regions where the nodes seem superfluous. The overall reduction of the number of nodes 
gives an optimal relation between precision and calculation cost (CPU time and memory constraint), 
by tracking the strong physical gradients within the flow field by higher grid point concentrations. 
This is especially important for unsteady flows. 

In this paper, inviscid compressible flow calculations from the ICASE Workshop on adapted grids 
are performed using dynamically auto-adaptive finite element triangular meshes. The mesh optimi- 
sation algorithms as well as the adaptation procedures are completely non-hierarchical, which allows 
more freedom for imposing optimisation strategies for obtaining admissible, regular grids, which can 
be geometrical, or structural and physical. The algorithm of the dynamical refinement /derefinement 
procedure is based on a certain number of basic algorithmic principles taking into consideration the 
particularities of local mesh refinement for finite element type generated meshes. A new anti-data 
structure has been adopted, where the successive subdivisions are performed independently of the 
former operations. 


LOCAL MESH REFINEMENT 


As the goal of mesh adaptation is to increase the accuracy of the solution process by locally 
enforcing the /i- adaptivity by usage of smaller discretisation cells, the first step in a mesh adaptation 
algorithm is to locally refine the mesh by adding new nodes according to some criteria, thus dimin- 
ishing local size of the concerned elements. The criteria used should be as close as possible to the 
error estimations of the underlying discretisation scheme. The principle of adaptive meshing is to 
uniformly equidistribute the error enhancing the overall convergence. 

Error estimates coming from the theory of finite-element simplex type meshes are based upon 
either a priori error estimates coming the governing equations (see e.g. [3]); those based on an a 
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posteriori error estimate where the computed residual of the solution R(u%) is used to define the error 
[4, 6]; and those which evaluate a combination of the derivatives of the computed variables (e.g. [7]). 
There is in fact a close relationship between such methods, based on physical gradients, and those of 
the second kind, [6], at least for hyperbolic transport equations. For the mesh adaptation procedures 
developed here, the latter strategy has proven to be robust and precise for tracking discontinuities for 
problems presenting strong and weak shocks, contact discontinuities and so on. For steady state flow 
resolution, local error criteria are based upon combinations of the L 2 average of V^Mach, VhDensity 
and S/hEntropy, with respect to a certain tolerance; whereas for unsteady flows, these are augmented 
by combinations of AhMach, AhDensity and AhEntropy. 

We describe here the techniques of mesh adaptation employed in the algorithm, further details 
can be found in [8, 9]. 

Symmetrical subdivision of triangles is obtained by adding a node in the middle of each side, 
rather than in the centre of each triangle, as only this method increases the precision of the numerical 
method, Figure 1. Indeed, from finite element theory, the interpolation is optimal when the triangle 
is as equilateral as possible within a certain metric. 



Figure 1: Symmetrical triangle subdivision. 


The local mesh refinement algorithm consists of testing the segments of an existing mesh whether 
they are to be split or not, which means that either 1, 2 or 3 new nodes will be created per cor- 
responding marked triangle. For geometrical considerations of maintaining non-acute angles, the 
creation of 2 new nodes is replaced by a symmetrical subdivision (Figure 1-right). The procedure 
performs refinement of triangles by subdividing them into 4 or 2 new elements. Division by two leads 
to zones of elements which may have a non-optimal distribution of the number of neighbours, leading 
to irregularities within the mesh. The close link between admissibility, regularity and optimality of 
an unstructured mesh to inherent properties of node neighbour numbers, has been fully exploited 
in the algorithms developed here, producing new constraints which are often more rigorous and less 
ad-hoc than existing ones where the adaptation criterion and its tolerance are the only reference 
points. For a 2D triangular mesh the optimal number of neighbours is 6. Typically, the two following 
refinement configurations should be avoided - Figure 2, and next section. 



Figure 2: Creation of nodes, with a non optimal node neighbour number. 

Their elimination is obtained by requiring that the the number of nodes that can be added within 
the 3 neighbouring elements should be greater than 2 and smaller than 5. This needs to be completed 
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by a geometrical criterion on the length relation between the smallest edge and the others in order 
to avoid creating highly stretched cells, (see next section). 

GEOMETRICAL CRITERIA OPTIMISATION 


After a refinement sweep, the new grid still “remembers” the old grid structure. To minimise such 
dependencies it is necessary to smooth globally the mesh as shown in Figure 3, using for instance 
the standard spring analogy : 

E a i x i 


Xi 


jefc(t) 


E % 


ctj — 1 and k(i) denotes the nodes neighbouring node i. The resulting mesh shows two specific 
behaviours : nodes having more than 6 neighbours tend to repulse their neighbours and those having 
less than 6 tend to attract them. 



Initial mesh 



Refined mesh 



Smoothed mesh 


Figure 3: Geometrical smoothing procedure. 


In order to reduce the above magnetic effect, a weight function related to the number of neighbours 
is introduced : 


aj — max [ 1 , 6 + (3 ( Hj — 6 ) ] 
with 


0 < /3 < 4 



Figure 4: Smoothing with a weight function. 


where Afj denotes the number of neighbours of node j. The result is illustrated in Figure 4. 

This optimisation could be more effective if nodes having more than 7 or less than 5 neighbours 
could be eliminated from the mesh. Such a requirement can be performed, by using the procedure 
of diagonal swapping with a neighbour node number minimising criterion. Let us denote by A/i, A/-J 
the number of neighbours for the vertices of a segment and A/ 3 , A /4 these numbers after swapping the 
segment as shown in Figure 5. Diagonal swapping is accepted if either : 
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M 3 + -A/4 < M\ -f- M 2 


or 


’ m 3 + m 4 

max (M3, M 4 ) 


= Mi + M 2 

< max (Mi, M 2 ) 




Figure 5: Diagonal swapping. 


Then the procedure allows to limit, in most cases, the maximum number of neighbours to 7 and 
thus increase the number of nodes having an optimal number of neighbours. 

To set the minimum number of neighbouring nodes, a technique of suppressing undesirable tri- 
angles can be used, which acts as shown in Figure 6. 



Figure 6: Cancellation of degenerate configurations. 


The algorithm used consists in marking segments having a vertex with less than 5 neighbours or 
with exactly 5 on each side. In both cases, the specific segment and their associate triangles, the 
shaded elements of Figure 6, can be suppressed. 

For the above square mesh, these techniques allow to regulate the optimal node neighbour number, 
and, combined with the weighted smoothing function, a considerable improvement in mesh regularity 
is achieved. 



Figure 7: Geometrical optimisation procedure. 


PHYSICAL CRITERIA OPTIMISATION 


The techniques of structural changes via moving nodes and diagonal swapping can also be applied 
depending on physical quantities. An improvement of the accuracy of the capture of a discontinuity 
can be obtained by aligning the edges. In Figure 8 the orientation of the edges were originally 
normal to a shock, which produces, on the left hand side, a relative thick shock. On the right hand 
side, a procedure based on diagonal swapping has been used, minimising the angles between the 
discontinuity and edge. 


222 





Standard refinement Refinement with edge alignment 

Figure 8: Improvement by aligning edges in the direction of the discontinuity. 


A modified spring analogy using a weight function related to a physical quantity can also be 
applied to obtain a squeezing property in regions of strong gradients. The weighted spring analogy 
is increased by the projection of neighbouring nodes onto the direction of the local gradient, and by 
weighting this projection by the local physical difference. This part of the algorithm corresponds to: 


Xi = 


5^ °b' d~ flij [ x i 4* 


nvftii 


a j + fo] 



Projection of node neighbours 
onto a gradient’s direction. 


where a j denotes the node number weight function, (3{j the local physical difference and V/?*/ || V#|| 
the gradient direction of the physical quantity. Here the density has been chosen as the specific 
physical quantity. 



Standard 




Refinement with squeezing 


Figure 9: Improvement by local squeezing around discontinuities. 


MESH DEREFINEMENT OR COARSENING 


In order to optimise the ratio precision and accuracy versus minimal number of discretisation 
nodes, dynamically auto-adaptive remeshing via local refinement and local derefinement (or coars- 
ening) is performed. Most local mesh derefinement techniques are based on a hierarchical data 
structure, which allows forward and backward scanning through the mesh, and maintains a history 
of previously added nodes. We have developed a new anti-hierarchical data structure; this algorithm 
enables nodes of an initial mesh to be removed, and also allows the use of the above structural 
optimisation techniques, which are incompatible with conventional hierarchical data structures. 
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This new algorithm is constructed as follows. The graph of points to be derefined is not straight- 
forward, so the inverse problem is solved. A list of nodes, called fixed points, which are never to be 
derefined is identified. These are either the singular nodes defining, for instance, the corners of the 
domain, symmetry points. .etc., or the nodes belonging to segments to be refined, or nodes having 
more than 6 neighbours. Then there are a number of possible fixed nodes allowing a maximal coars- 
ening which come from either boundary or internal configurations. To find these, first the border of 
the domain is scanned so that each second node becomes a fixed point, and finally the configuration 
of Figure 10 is searched through the mesh. Thus all triangles containing no fixed point have at least 
three fixed neighbouring nodes, as shown in Figure 10. 



Figure 10: An element defined with 3 removable nodes should be surrounded by 3 fixed nodes. 


Thus the grid becomes coloured according to the nodes to be removed and those to be kept in 
the mesh. The derefinement can then be done by performing successively the operation presented in 
Figure 11. It is necessary to treat firstly the element with 3 nodes to be invalidated, then the cases 
with 2 nodes to be removed, and finally the case having only one node which disappear. The dark 
shaded elements of Figure 11 vanish during the procedure and the light shaded triangles change their 
shape. 



If ^ JJ. 



Remove triangles with Remove triangles with Remove triangles with 

3 marked nodes. 2 marked nodes. 1 marked node. 

Figure IT. Derefinement stages. 
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ALGORITHM 


The whole algorithm becomes thus : first, a list of nodes to be invalidated is defined, followed by 
their destruction, even if they belong to the initial macro-mesh. Then the phase of local enrichment 
(adding new nodes) is performed by adding them in the middle of segments. This is followed by both 
structural optimisation procedures, (first the physical diagonal swapping and then the geometrical 
ones), and finally the mesh is smoothed with the weighted spring analogy to obtain a regular, confor- 
mal and admissible mesh. For steady state simulations the squeezing procedure is performed after a 
specified convergence on the last two meshes. The total number of different meshes created depends 
upon the physical solution, for steady state calculations, between 4 to 8 different stages of remeshing 
are performed. The whole procedure is directly integrated into the flow solver giving a fully dynamic 
auto-adaptive mesh algorithm. For unsteady calculations, the remeshing depends greatly upon the 
speed of the transient states, and for a calculation cycle, there can be as many as 400 or more new 
meshes created. The increase of CPU time introduced by these procedures is the order of 2 to 5 
% for steady state calculations, and is inferior to 30 % for unsteady ones. The internal accounting 
of memory requirements are managed by a dynamical memory manager and a archival data base 
structure; an updated mesh replaces the former one which requires only a minimal extra memory 
allocation. 


APPLICATIONS FOR EXTERNAL FLOWS OVER PROFILES 

The flow solver used for all these applications is based upon an equivalent Galerkin finite volume 
approximation on the dual control volumes of the PI triangular simplex. A Jameson type centred 
scheme is employed with artificial dissipation, [9, 8]. 

Transonic flow over a NACA0012 at Moo = 0.80, a = 1.25° 

This is the AGARD 01 test case [1], concerning transonic flow over a NACA profile. Its particular- 
ities are the presence of a weak shock on the windward side which is a good challenge for testing shock 
capturing criteria for local refinement and squeezing in this region, where the gradients are not as 
strong as those that can be observed, for instance, with another very standard test case, Moo = 0.85 
a — l.°, (AGARD 02). The shock on the leeward side is a much sharper discontinuity, and thus, is 
easier to localise. Since the test case concerns a non-zero angle of attack flow over a profile in a 
bounded domain, a non-zero circulation is induced on the profile, it is necessary to correct the outer 
boundary condition in order to simulate as well as possible an infinite domain [10]. The sensibility of 
the values of the lift coefficient to this induced circulation effect is high. The outer boundary should 
be chosen as far as possible, and be of a suitable dimension, as well as a “vorticity” correction which 
is applied to the infinite boundary conditions. For transonic flows, the infinite boundary values cor- 
respond to three entering characteristics and one outgoing one. The outgoing component can allow 
the generated circulation to influence the infinite boundary values of the solution. 

Several different geometrical shapes for the domain boundaries were tested, as well as varying 
profile to infinite boundary lengths, as a function of the chord length. For a fixed outer boundary 
distance, the best results were obtained by a circular outer boundary. The aerodynamical coefficients 
are tabulated below, where a significant variation is found even by taking the precaution of applying 
the vorticity correction. The values stabilize after as distance of 1000 chord lengths, (up to 10000 
chord lengths were tested). They correspond well to the references values established by AGARD. 

The solution adaptive mesh procedures presented in this paper were applied, and despite the 
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Table 1: Flight coefficients for NACA0012 at M m = 0.80, a = 1.25° 


Outer boundary 

Ci 

c d 

c m 

X U pp er 

X lower 

10 chords 

.3375 

.0224 

-.0351 

.6274 

.3441 

30 chords 

.3495 

.0232 

-.0371 

.6327 

.3422 

100 chords 

.3549 

.0235 

-.0382 

.6353 

.3382 

1000 chords 

.3569 

.0237 

-.0387 

.6363 

.3357 


difficulties of capturing the weak windward side shock, a highly satisfactory adaptive mesh was 
obtained. For this test case and the next one, where shocks and expansion fans are present within 
the flow field, a combination of gradients and differences of physical variables were chosen to provide 
the adaptation criteria as explained in the previous sections. All six distinct meshes were generated 
during the process. The initial mesh was of 1704 nodes and 3332 elements, the final one of 8136 
nodes and 16078 elements, Figures 12. 






Figure 12: Evolution of the meshes and Iso C p lines for a NACA0012 at = 0.80, and a = 1.25°. 

A partial view is given below, Figure 13, showing clearly the various types of adaptation per- 
formed, local refinement, regularisation, alignment, squeezing, and the corresponding pressure coef- 
ficient body profiles are presented to demonstrate the precision of solver on this adapted mesh. 
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Figure 13: Mesh Details, Iso C p lines (AC P = 0.1) and C p body coefficients. 


Transonic flow over a NACA0012 at Moo = 0.95, a = 0° 

The second test case over the NACA 0012 presented here corresponds to the AGARD 03 test 
case. Despite the zero angle of attack for this high transonic flow, there is again a great sensibility of 
the solution to the position of the outer boundary. Indeed, the test case gives rise to a fish-tail shock 
structure, with an oblique shock attached to the trailing edge and a normal shock emanating from 
the intersection of the sonic line with the trailing edge oblique shock (see Figures 14). The distance 



1704 nodes, 3332 elements 39586 nodes, 78864 elements 54472 nodes, 108628 elements. 



Figure 14: Evolution of meshes and Iso C p lines, for a NACA0012 at M TC = 0.95, a = 0.0° 
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Figure 15: Meshes and C p solutions for different outer boundary distances (10, 100, 10000 chords). 


of this normal shock to the trailing edge is extremely sensitive not only to the position of the outer 
boundary, but also to the degree of adaptation around the body itself due to the strong expansion 
waves which are created along the profile. 

The outer boundary was placed successfully at 10, 30, 100, 1000 and 10000 chord lengths from the 
body. Whereas the flight coefficients converged after 100 chords, the normal shock position continues 
to change throughout the series, increasing with increasing outer boundary distance. The degree of 
adaptation for these comparisons (Figures 15) was kept constant at a level equivalent to those of the 
final mesh shown in the Figures 14. The evolution of the refinements for the case with 100 chords 
are given in the Figures 14. 

The convergence of the flight coefficients and the distance Xs of the normal shock to the trailing 
edge (normalised to the chord length) are tabulated below. As mentioned above, the normal shock 
distance evolves with the outer boundary distance, as has also been observed in the literature, [1]. 
The orders of magnitude of the drag coefficient are in agreement with the references, however the 

Table 2: Flight coefficients for NACA0012 at = 0.95, a = 0.0° 


Outer boundary 

C, 

c d 

Cm 


10 chords 

.0007 

.1090 

-.0002 

1.280 

30 chords 

-.0003 

.1091 

.0002 

2.515 

100 chords 

.0004 

.1092 

-.0001 

3.125 

300 chords 

-.0009 

.1091 

.0001 

3.179 

1000 chords 

.0003 

.1091 

.0000 

3.195 

10000 chords 

-.0004 

.1091 

.0000 

3.231 
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normal shock stand-off distance is found to be closer than the references for a low number of chord 
distances for the outer boundary. The calculations presented within this monograph however were 
often on non-adaptive grids, and for a relatively low number of chord lengths for the outer boundary 
distance. 


Subsonic flow over a multi element airfoil at M m = 0.2, a = 0.0° 

This test case, corresponds to a four-element airfoil of Suddhoo and Hall, obtained by applying the 
Karman-Trefftz mapping function, [2]. The flow conditions are subsonic with zero angle of attack. 
The adaptation criteria here were taken to be based upon simply the gradient of the local Mach 
number, as there are no shocks present within the flow field. Again, this test case is a steady state 
calculation. Despite the absence of discontinuities the mesh adaptation algorithms provided very 
regular final meshes, starting from an initial coarse and non-optimal mesh, Figures 16. 



Figure 16: Evolution of the meshes and Iso C p lines, (A C p = 0.3.) for the Multi-Element airfoil test 
case, for a subsonic flow of = 0.2, a = 0.0° 


Airfoils presenting Non uniqueness of the solutions for the Euler equations 

In an AIAA paper concerning airfoils that can occur during an optimum design procedure, where 
the surface splines are perturbed, A. Jameson found the possibility of obtaining two distinct solution 
branches of the Euler equations, whereas the meshes used were extremely fine, and convergence was 
pushed to a maximum, in order to eliminate possibilities of anomalies due to a non-entropy preserving 
solution [5]. In fact, the Euler equations only admit weak solutions, and the only admissible ones 
must preserve the mathematical entropy condition. The airfoils studied here show a hysteresis effect 
in the lift /incidence polar when varying the angle of attack back and forth. The solutions generated 
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have two separate supersonic zones and shock waves for a particular angle of attack, and then by 
increasing or decreasing the lift angle around some critical value, a single supersonic zone is generated 
with a shift in the lift coefficient. It is extremely important to start this varying a process with an 
extremely fine and regular grid around the airfoil. The freestream Mach number is fixed at .78, 
and an initial series of computations is performed upto complete convergence for increasing angle of 
attack from zero to —.70, (up to 30000 time step iterations or more). The different lift coefficients 
are noted. Then for the solutions around a Cl of .6, angle of attack around —.43, a second series 
are initialised by taking the solution of an angle of- attack slightly inferior. As the mesh adaptation 
is dynamically linked to the solver, the essential characteristics of the changing solution are taken 
into account, and each calculation is thus made on its own Specific adapted mesh. Several different 
stages of resolution are illustrated in the Figures 19. The hysteresis effect is plotted in Figure 18. 
The values do not correspond exactly to those of Jameson, as the initial definition of the airfoil was 
not sufficiently detailed to obtain an equivalent initial shape. The transient stages are well captured 
by the dynamical mesh method, and the non-uniqueness can be established within a certain margin 

.75 
.70 
.65 
Q .60 
.55 
.SO 
.45 

-.70 -.65 .60 -.55 -SO -.45 -.40 -.35 -.30 -.25 -.20 

aipha 

Figure 18: Lift versus angle of attack polar for Mach^ = .78 
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Figure 19: Different Meshes and Mach isolines/body plots obtained for M <*> — .78, a — —.47°, 
corresponding to hysterisis effect of the lift polar. 


of Mach, angle of attack and lift. 


CONCLUSIONS 

A new “anti 1 - data structure unstructured mesh adaptation procedure has been presented enabling 
the use of structural optimisation, together with freedom for derefinement of any level of mesh 
transformation. The improvement, for steady Euler flows solutions, in number of nodes and CPU 
time requirement can be up to a factor of 30, and the algorithms developed have also proved to 
produce very accurate results for unsteady cases. The methods have been applied to a number of 
test cases, and have proved to be a valuable design tool for aerodynamic investigation. The solver 
and the dynamical mesh algorithm are intimately coupled, their implementation on a distributed 
computational platform was undertaken in a Master-Slave environment between a CRAY YMP and a 
CRAY-T3D. The efficiency of such a technique is extremely promising, and in vue of the minimisation 
of the communication phases by performing the adaptation and the partitioning on the Master, proves 
to be an interesting alternative to a complete MPP algorithm of adaptation and solution process, 
which necessitates sorting graphs, redistribution searches and high communication patterns. 
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